Hydrogen bonds and van der Waals forces in ice at ambient and high pressures 
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The first principles approaches, density functional theory (DFT) and quantum Monte Carlo, have 
been used to examine the balance between van der Waals (vdW) forces and hydrogen (H) bonding 
in ambient and high pressure phases of ice. At higher pressure, the contribution to the lattice energy 
from vdW increases and that from H bonding decreases, leading vdW to have a substantial effect 
on the transition pressures between the crystalline ice phases. An important consequence, likely 
to be of relevance to molecular crystals in general, is that transition pressures obtained from DFT 
functionals which neglect vdW forces are greatly overestimated. 



Water-ice, the most common molecular solid in nature, 
exhibits a rich and complex phase diagram. At present 
this includes ice Ih, ice Ic and 14 other crystalline ice 
phases [2 |2] . Although the phase diagram of ice up to 
around 2 GPa is well-established experimentally, under- 
standing of the subtle balance of intermolecular interac- 
tions which give rise to this richness is incomplete. In par- 
ticular, the relative contribution of hydrogen (H) bonding 
and van der Waals (vdW) dispersion forces to the cohe- 
sive properties of the various crystalline ice phases is still 
not understood. From a theoretical perspective, the var- 
ied densities of experimentally characterized ice crystal 
structures, with the molecules fixed on well-defined lat- 
tices, affords an excellent opportunity to quantify and 
assess the nature of these intermolecular interactions. 

Computer simulation techniques have proved instru- 
mental in understanding ice (e.g. [3l-[T4]). In particular, 
density-functional theory (DFT) with generalized gradi- 
ent approximation (GGA) functionals has been widely 
applied. Certain GGAs describe the ambient pressure 
ice Ih phase reasonably well [9] and predict the proton 
order-disorder phase transition temperatures between ice 
Ih and XI and ice VII and VIII in good agreement with 
experiments [B]. However, it is known that GGAs suffer 
deficiencies when vdW forces are important and indeed it 
has been suggested that this is the reason certain GGAs 
produce a liquid water density about 15-20% less than 
experiment |15H17) or fail to predict the correct ground 
state structure for some small water clusters [Hj . Despite 
the considerable body of DFT work on ice phases, the 
role of vdW forces and H bonding has not been system- 
atically examined in anything other than ice Ih [191 EQ] . 
However, with recent developments (e.g. pTk23j ). it is 
now possible to tackle this issue head on and estimate the 
importance of vdW forces and H bonding in the various 
phases of ice. 
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TABLE I. Absolute lattice energies (omitting zero point en- 
ergy effects) of ice Ih, II, and VIII. Relative energies com- 
pared to ice Ih (AUo) are given in parenthesis. All values are 
in meV/HaO. 





Ih 


II 


VIII 


Exptj^] 


-610 


-609 (1) 


-577 (33) 


DMC 


-605 ±5 


-609 ±5 (-4) 


-575 ±5 (30) 


PBE 


-636 


-567 (69) 


-459 (177) 


PBEO 


-598 


-543 (55) 


-450 (148) 


PBEO-HvdW'^^ 


-672 


-666 (6) 


-596 (76) 



^ Ref. | 24| . with zero point energy contributions removed. 



Here we report an extensive series of first principles 
studies aimed at better understanding the role of vdW 
and H bonding in ice. This includes DFT calculations 
with and without a treatment of vdW forces on the am- 
bient pressure phase of ice, ice Ih, and all the proton 
ordered phases, namely, in order of increasing pressure, 
ice IX, II, XIII, XIV, XV, VIII. We also report diffusion 
quantum Monte Carlo (DMC) calculations for ice Ih, II, 
and VIII as reference values to complement the experi- 
mental results. DMC is highly accurate for weak inter- 
actions including vdW bonded systems and represents 
the state-of-the-art for electronic structure simulations 
of solids. From this work we find that the contribution 
of vdW to the cohesive properties of ice increases as one 
moves to the higher density phases, whereas the contri- 
bution from H bonds decreases. As a result, vdW plays 
a crucial role in determining the relative stabilities and 
phase transition pressures in ice. The results presented 
here are likely to be of relevance to understanding in- 
termolecular interactions in water in all its condensed 
phases as well as to structural searches for (novel) high 
pressure ices and to molecular crystals in general. 

We start by discussing lattice energies, which are ob- 
tained by subtracting the total energy of the ice unit 
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FIG. 1. Relative lattice energies, AUo (a), and volumes, AVo (b), of the high pressure ice phases with respect to the lattice 
energy of ice Ih obtained with the methods indicated, (c) Transition pressures (Ptr) from ice Ih to the various high pressure 
phases. 



cell containing n H2O molecules from the total energy of 
n isolated H2O molecules. In this context Whalley's ex- 
trapolations of the experimental finite temperature and 
pressure phase coexistence lines to zero temperature and 
pressure are extremely valuable (23]. For the proton or- 
dered phases Whalley considered, these indicate that ice 
IX, II, and VIII are less stable than ice Ih by only 3.5 ± 
0.8, 0.6 ± 1.0, and 33 meV/HaO, respectively [24l. These 
values agree well with DMC. Specifically, DMC predicts 
that ice VIII is 30 ± 7 meV/H20 less stable than ice Ih. 
Similarly, the near energetic degeneracy between ice Ih 
and ice II is captured with DMC, and within the DMC 
error bars, ice Ih and ice II are equally stable. From Ta- 
ble I it can also be seen that DMC lattice energies are 
in very good agreement with experiment as well. Overall 
this gives us confidence in the quality of Whalley's ex- 
perimental extrapolations and the accuracy of the DMC 
calculations. This therefore provides an excellent basis 
for exploring the role of H bonds and vdW in ice with 
DFT. 

We now discuss the results obtained with FEE [55], 
one of the most widely used functionals. For ice Ih, PBE 
yields a reasonable lattice energy of about 640 meV/Il2 0, 
an overestimate of around 30 meV/Il20 which is con- 
sistent with previous work However, Table I and 
Fig. 1(a) reveals a severe deterioration in the perfor- 
mance of PBE for the higher density phases. Whereas 
experiment and DMC suggest that the energy difference 
between ice Ih and the least stable ice VIII phase is 
about 30 meV/H20, PBE gives ~180 meV/H20 differ- 
ence. This is mainly because PBE underestimates the 
stability of the high pressure phases, suggesting that at- 
tractive interactions, more important in the higher den- 
sity ice phases, are not captured accurately with PBE. In 
addition to substantially overestimating the energy dif- 
ferences between the various phases, the changes in vol- 
umes (AVo) upon going from ice Ih to the high pressure 
phases are also underestimated with PBE (Fig. 1(b)) 
and critically the transition pressures (Ptr) — obtained 
from Ptr = — AL/q/AVo — are ~5-15 times larger than 



experiment (Fig. 1(c)). 

Seeking to understand why PBE performs so poorly for 
the high pressure ice phases, we have considered various 
potential sources of errors. First we analyzed the dipole 
moment and polarizability of an isolated water molecule. 
PBE, like other GGAs, overestimates the polarizability 
of an isolated water molecule by ~10% compared to ex- 
periment, which is related to any overly delocalized elec- 
tron density and a too small HOMO-LUMO gap. In- 
corporation of a fraction of Hartree-Fock exchange in 
to the GGAs is an established approach for alleviating 
this problem. In particular, a popular PBE-based hybrid 
functional, PBEO g^, widens the HOMO-LUMO gap of 
the isolated water molecule by ^-^40% and provides more 
accurate polarizabilities and interaction energies within 
a variety of water clusters [55^50] . When applied to ice, 
the lattice energies obtained with PBEO (Table I) arc in- 
deed improved over PBE, however, the energy difference 
between ice Ih and ice VIII is still about four times larger 
than the experimental value. In addition, as with PBE, 
upon going from ice Ih to the high pressure phases the 
transition pressures are about an order of magnitude too 
large compared to experiment (Fig. 1(c)). 

Turning to the role of H bonding, the relative H 
bond strength has been estimated from shifts in the 0-H 
stretching frequencies, which is a simple and widely used 
measure of H bond strength (e.g. [21]). Specifically, the 
softening (red-shift) of the intramolecular 0-H stretch- 
ing frequencies in the various phases is compared to the 
(average) 0-H stretching frequency of an isolated water 
monomer [32' . Based on the available experimental fre- 
quencies for several of the phases ^33.-^55^, H bonds get 
weaker with increasing pressure. This established, but 
not necessarily obvious, result arises because the nearest 
neighbor water-water distances get larger as one moves 
from ice Ih with its open ring structure to the more com- 
plex higher pressure phases. In line with previous cal- 
culations [36], PBE reproduces the trend in frequency 
changes with pressure but the red-shift (i.e., strength) 
compared to the the isolated water molecule is signifi- 
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cantly overestimated compared to experiment, particu- 
larly for ice Ih (see Fig. 2(a)). PBEO predicts frequency 
shifts in much better agreement with experiment (Fig. 
2(a)). However, as we have seen, the relative stabilities 
of the ice phases with PBEO are significantly in error, 
which clearly suggests that H bonds are not the only im- 
portant interaction in the high density phases. 

We now consider the influence of vdW interactions 
with the scheme of Tkatchenko and Scheffler (referred to 
as vdW"""^), in which an additional Cg/i?^ tail is added to 
the DFT total energy, with the Cg coefficients calculated 
as functionals of the electron density [22] . Dramatic im- 
provements in the relative energies of the various phases 
and phase transition pressures are observed (Fig. 1). In 
particular, upon going from PBEO to PBEO+vdW''"^ the 
energy difference between ice Ih and ice II is reduced to 
just 6 meV/Il20 and, likewise, ice VIII is now only 76 
meV/Il20 less stable than ice Ih. As a consequence, with 
vdW the ice Ih to ice II transition pressure is reduced by 
more than one order of magnitude and is in good agree- 
ment with experiment. Similarly, the ice Ih to ice VIII 
transition comes within a factor of two of experiment. 
The substantially improved transition properties when 
vdW is accounted for results from a strong dependence 
of vdW on density, as shown in Fig. 2(b). Clearly as 
the density of the ice phases increases, so too does the 
vdW contribution to the lattice energy and indeed for the 
highest density phase (ice VIII) the vdW contribution is 
about twice what it is in ice Ih. Overall the decrease of 
H bond strength and the increase of vdW in the high 
pressure phases clearly shows a significant enhancement 
of relative contribution of vdW over H bonding interac- 
tions for the cohesive properties of the high pressure ice 
phases. 

Our findings regarding the importance of vdW for the 
high-pressure phases of ice also hold when using function- 
als based on the vdW-DF approach of Dion et al. pTll37]. 
The phase transition pressures obtained with such func- 
tionals are similar to those obtained with the TS scheme. 
Both approaches neglect the non-additive many-body 
vdW energy beyond the pairwise approximation. We 
have found that the many-body vdW energy plays a mi- 
nor role for the different phases of ice by using an exten- 
sion of the TS scheme [38j- Also we find that quantum 
effects based on zero point energies play a minor role in 
determining the relative stabilities of the various phases. 
Specifically, the zero point energies of the different phases 
(calculated with PBE and the harmonic approximation) 
differ by <10 meV/HsO. 

We now discuss why vdW interactions play a crucial 
role in determining the relative stabilities of the various 
ice phases. For any two non-bonded atoms vdW forces 
result from an induced dipole-induced dipole interaction 
whose leading term varies as CqR^^ . At a given i?, vdW 
will increase if the polarizability of the atoms increases, 
i.e., Cg becomes larger and/or the number density of the 
atoms increases. Our calculations with the TS approach 
show that the Cg coefficients do indeed increase upon go- 




Density [g/cm^] Density [g/cm'] 

FIG. 2. (a) H bond strength index [32] and (b) vdW energy 
contributions plotted as a function of the experimental den- 
sities of the ice phases at zero pressure. Experimental values 
are taken from [33ri35) . 



ing from ice Ih to ice VIII, by 24% for H and 6% for O. 
However, the major effect that leads to an increase in the 
vdW interactions in the high density phases is simply the 
higher packing of water molecules. This can be seen by 
comparing the 0-0 radial distribution functions of ice 
Ih with e.g. ice II and ice VIII (Fig. 3(a)). Although 
ice Ih possesses the shortest nearest neighbor 0-0 dis- 
tances, the structure is open and there is a large gap of 
~2 A between the first and second coordination shells. 
In contrast, in the high density phases the second and 
subsequent coordination shells appear at much shorter 
0-0 separations. In fact in ice VIII, which is comprised 
of two interpenetrating sub-lattices, the first and second 
coordination shells fall almost on top of each other (with 
the shortest 0-0 distances associated with non H bonded 
contacts). The higher packing, particularly in the ca. 3 to 
6 A regime is reflected by the integrated number of neigh- 
bors versus 0-0 distance shown in Fig. 3(b)). Overall 
these additional molecules at short (ice VIII) and inter- 
mediate (ice II) distances lead to enhanced vdW interac- 
tions in the high pressure phases. Fig. 3(c) provides a 
more quantitative basis for the above argument, by show- 
ing the total vdW interactions — as obtained from the 
TS scheme — for ice Ih, II and VIII in the range of ca. 3 A 
to 6 A. It can be seen that beyond about ca. 3 A there are 
more substantial contributions from vdW from ice II and 
VIII than ice Ih. Although vdW is generally considered 
to be a long range interaction, the dominance of vdW in 
this ca. 3 A to 6 A regime is to be expected given that 
the sum of the vdW radii of O and H in the condensed 
phase (calculated with TS scheme) is only <3 A. Note 
that in fact the difference in the vdW contribution be- 
tween the various phases is essentially converged to the 
periodic limit after about 8 A. 

In conclusion, we have performed an extensive first 
principles study of ice at ambient and high pressures. 
From the ambient to the high pressure phases the con- 
tribution to the lattice energy arising from vdW forces 
monotonously increases. This has a signiflcant impact 
on the phase transition pressures, as exemplifled by cal- 
culations with conventional xc functionals where vdW 
forces are not accounted for, and where transition pres- 
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FIG. 3. (a) O-O radial distribution functions {g{r)) of ice Ih, 
ice II, and ice VIII. (b) Integrated number of neighbors {N{r)) 
and (c) vdW contributions (obtained from the TS scheme) as 
a function of O-O distance for the same three phases. 



sures more than one order of magnitude larger than ex- 
periment are obtained. By accounting for vdW forces, 
the phase transition pressures more closely agree with 
experiment. This finding provides a new physical insight 
and is also relevant to DFT-based structural searches for 
novel ice polymorphs |14|, I39j , by implying that a reex- 
amination of the high pressure region of the ice phase 
diagram with a vdW corrected DFT approach would be 
worthwhile. Whilst the focus of the current study has 
been on understanding the role of H bonding and vdW 
forces in ice, rather than on improving the current state- 
of-the-art in DFT simulations for water, it is nonetheless 



worthwhile to stress that there certainly remains scope 
for further improvements within DFT in terms of transi- 
tion pressures, absolute lattice energies and volumes for 
ice. In particular, there is scope for improving the lat- 
tice energy of ice VIII obtained with vdW"'"'^, with the 
current underestimation in the lattice energy possibly as- 
sociated with the radically different structure of ice VIII 
compared to the other phases considered. The strong 
dependence of vdW interactions on density suggests that 
even at ambient pressure, vdW forces will also be impor- 
tant to liquid water which has a ^8% higher density than 
ice Ih. Indeed this is consistent with a number of recent 
DFT studies [Tni[T7]. Although vdW is often associated 
with so-called "sparse" matter we have shown that it at- 
tains greater significance at high density. This somewhat 
counterintuitive result arises because vdW forces mainly 
enhance interactions between molecules at medium range 
(i.e., second and third coordination shells at 3 to 6 A 
separations). Therefore, analogous effects are expected 
for water in other environments such as in confined ge- 
ometries, interfaces, and clathrates, and indeed for other 
hydrogen bonded molecular crystals at high pressure. 
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